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ABSTRACT 

We report 21-year timing of one of the most precise pulsars: PSR J1713+0747. Its pulse times of arrival 
are well modeled by a comprehensive pulsar binary model including its three-dimensional orbit and a noise 
model that incorporates short- and long-timescale correlated noise such as jitter and red noise. Its timing 
residuals have weighted root mean square ~ 92 ns. The new data set allows us to update and improve previous 
measurements of the system properties, including the masses of the neutron star (1.31 ±0.11 M 0 ) and the 
companion white dwarf (0.286 ±0.012 Mg) as well as their parallax distance 1.15 ± 0.03 kpc. We measured 
the intrinsic change in orbital period, f£ nt , is -0.20 ±0.17 ps s _1 , which is not distinguishable from zero. 

This result, combined with the measured P^' of other pulsars, can place a generic limit on potential changes 
in the gravitational constant G. We found that G/G is consistent with zero [(-0.6± 1.1) x 10 ~ 12 yr _1 , 95% 
confidence] and changes at least a factor of 31 (99.7% confidence) more slowly than the average expansion 
rate of the Universe. This is the best G/G limit from pulsar binary systems. The P/" of pulsar binaries can 
also place limits on the putative coupling constant for dipole gravitational radiation kd = (—0.9 ± 3.3) x 10 -4 
(95% confidence). Finally, the nearly circular orbit of this pulsar binary allows us to constrain statistically the 
strong-field post-Newtonian parameters A, which describes the violation of strong equivalence principle, and 
tt 3 , which describes a breaking of both Lorentz invariance in gravitation and conservation of momentum. We 
found, at 95% confidence, A < 0.01 and ±3 < 2 x 10 -20 based on PSR J1713+0747. 

Subject headings: pulsars: individual (PSR J1713+0747) — Radio: stars — stars: neutron — Binaries:general 
— gravitation - relativity 


1 . INTRODUCTION 

We present 21-year timing of the millisecond pulsar (MSP) 
J1713+0747, which was discovered in 1993 (Foster et al. 
1993). It is one of the brightest pulsars timed by the North 
American Nanohertz Observatory for Gravitational Waves 
(NANOGrav; McLaughlin 2013; Demorest et al. 2013), and 
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has the smallest timing residual of all NANOGrav pulsars 
(Demorest et al. 2013). The timing analysis reported in this 
paper will be incorporated into future pulsar timing array 
projects. Timing observations of this pulsar were reported 
previously in Camilo et al. 1994, Lommen & Backer 2001, 
van Straten & Bailes 2003, Splaver et al. 2005, Hotan et al. 
2006 and Verbiest et al. 2009. We report new results that arise 
from a significant extension of the timing baseline and use of 
wide-bandwidth, high resolution instruments that allow us to 
model and account for pulse time variations due to dispersion 
in the interstellar medium (ISM) (Section 3.2) and small dis¬ 
tortions of the pulsar’s magnetosphere (Section 3.4). 

MSPs are very stable rotators due to their enormous angu¬ 
lar momentum. PSR J1713+0747 is an MSP residing in a 
wide binary orbit with a white dwarf companion (Section 3). 
The pulse arrival times of the pulsar are well fit by a binary 
model with a nearly circular orbit. The masses of the binary 
components can be inferred through the measurement of the 
mass function and Shapiro delay (Splaver et al. 2005). The 
system’s distance is well-measured through a timing parallax. 
We detect a changing projected orbital semi-major axis due 
to the orbit’s proper motion on the sky. Through the rate of 
change of projected semi-major axis, we can infer the orien¬ 
tation of the orbit in the sky. This is one of the few bina¬ 
ries in which the 3D-orientation of the binary orbit can be 
completely solved. Using 21 years of data, we refine the pre¬ 
viously published measurements of these orbital parameters. 
We observed apparent variation of the binary orbital period 
due to the Shklovskii effect and Galactic differential acceler¬ 
ation. We also find a stringent constraint on the intrinsic vari¬ 
ation of the orbital period, which enable us to test alternative 
theories of gravitation. 

The stability and long orbital period of the 
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PSR J1713+0747 binary make it an excellent laboratory 
for observing the time variation of Newton’s gravitational 
“constant” G. This interesting conjecture of G varying on 
a cosmological timescale was first raised by Dirac (1937) 
based on his large-number hypothesis. Later this become a 
prediction of some alternative theories of gravitation. For 
example, the scalar-tensor theory (Jordan 1955, 1959; Fierz 
1956; Brans & Dicke 1961) modifies Einstein’s equation of 
gravitation by coupling mass with long-range scalar-tensor 
fields, and predicts that, as the universe expands, the scalar 
field will also change, causing the effective gravitational 
constant to vary on the cosmological timescale. Similar 
ideas were also revisited by gravitational theories involving 
extra dimensions (Marciano 1984; Wu & Wang 1986). 
PSR J1713+0747 is likely the best pulsar binary for testing 
the constancy of G thanks to its high timing precision and 
long orbital period. Using timing results reported in this 
paper, we found a stringent generic upper limit on G (see 
Section 4.3 for details). 

The PSR J1713+0747 binary is also an excellent laboratory 
for testing the strong equivalence principle (SEP) and the pre¬ 
ferred frame effect (PFE; in this work, we constrain the pu¬ 
tative post-Newtonian parameter 6:3 that characterizes a spe¬ 
cific type of PFE; for details see Section 4.4). The violation of 
SEP or the existence of non-zero 6:3 could lead to potentially 
observable effects which cannot be accounted in the context 
of GR, such as forced-polarization of the binary orbit. Some 
of these effects are discussed in Freire et al. (2012a) and Will 
(2014). We can put stringent generic constraints on alternative 
theories by observing low-eccentricity pulsar binary systems. 
Section 4.4 presents the constraint on the violation of SEP 
and the significance of <+3 from our 21 years of observation of 
PSR J1713+0747. 

2 . OBSERVATIONS 
2.1. Data acquisition systems 

Our data set consists of pulse timing observations of PSR 
J1713+0747 at the Arecibo Observatory, from 1992 through 
2013, and at the Robert C. Byrd Green Bank Telescope 
(GBT), from 2006 through 2013, using several generations of 
data acquisition systems as described below (Table 1). The 
first 12 yr of these data (1992 August through 2004 May) 
were previously reported by Splaver et al. (2005) using the 
Mark III, Mark IV, and Arecibo-Berkeley Pulsar Processor 
(ABPP) systems. In this work, we incorporate an extension 
of the Mark IV data set (reduced using the same process as 
Splaver et al. 2005); data collected at Arecibo from two newer 
systems; and data collected at the GBT. 

The earliest observations (Mark III) were made using a sin¬ 
gle receiver operating over a limited bandwidth and hence in¬ 
trinsic pulsar behavior cannot be separated from effects of the 
ISM. Subsequent observations at both telescopes were made 
using two receivers at widely spaced frequencies (1410 and 
2380 MHz at Arecibo; 800 and 1410 MHz at the GBT) to 
allow for separation of these effects. 

The earliest observations (1992-1994) used the Prince¬ 
ton Mark III (Stinebring et al. 1992), which collected dual¬ 
polarization data with a filter bank of 32 spectral channels 
each 1.25 MHz wide. Observations between 1998 and 2004 
used the Princeton Mark IV (Stairs et al. 2000) instrument 
and the ABPP (Backer et al. 1997) system in parallel. The 
Mark IV system collected 10-MHz passband data using 2-bit 
sampling. The data were coherently dedispersed and folded at 


the pulse period offline. The ABPP system sampled voltages 
with 2-bit resolution and filtered the passband into 32 spec¬ 
tral channels (1.75 MHz per channel and 56 MHz in total for 
1410-MHz band; 3.5 MHz per channel and 112 MHz in to¬ 
tal for 2380-MHz band), and applied coherent dedispersion to 
each channel using 3-bit coefficients. 

From 2004 to 2011/12, pulsar data were collected with the 
Astronomical Signal Processor (ASP; Demorest 2007) and 
its Green Bank counterpart GASP (Demorest 2007). The 
(G)ASP systems recorded 8 -bit sampled ~64 MHz band¬ 
width data and applyed real-time coherent dedispersion and 
pulse period folding. The resulting data contain 2048-bin full- 
Stokes pulse profiles integrated over 1-3 minutes. When ob¬ 
serving with ASP we used 16 channels each 4-MHz wide be¬ 
tween 1440 and 1360 MHz, and 16 channels between 2318 
and 2382 MHz. When observing with GASP we used 12 
channels between 1386 and 1434 MHz, and 16 channels be¬ 
tween 822 and 886 MHz. The J1713+0747 ASP/GASP data 
were also reported in Demorest et al. (2013). 

We started using the Green Bank Ultimate Pulsar Process¬ 
ing Instrument (GUPPI; DuPlain et al. 2008) for GBT obser¬ 
vations in 2010 and its clone the Puerto-Rican Ultimate Pulsar 
Processing Instrument (PUPPI) for Arecibo observations in 
2012. GUPPI and PUPPI use 8 -bit sampling in real-time co¬ 
herent dedispersion and pulse period folding mode and pro¬ 
duce 2048-bin full-Stokes pulse profiles integrated over 10 
second intervals. When observing with the 800-MHz receiver 
at the GBT, we use GUPPI to collect data from 62 spectral 
channels each 3.125 MHz wide, covering 724-918 MHz in 
frequency. With the L-band receiver, GUPPI uses 58 spectral 
channels each 12.5 MHz wide, covering the 1150-1880 MHz 
band. When observing with the L-band receiver at Arecibo, 
we use PUPPI to collect data from 1150-1765 MHz using 
50 spectral channels each 12.5 MHz wide. When observing 
with the S band, PUPPI takes data from 1770-1880 MHz and 
2050-2405 MHz using 38 spectral channels each 12.5 MHz 
wide. The spectral bandwidth and resolution provided by 
GUPPI and PUPPI are crucial for resolving the pulse profile 
evolution in frequency described in Section 3.4. 

2.2. Arrival time calculations 

We combine the pulse times of arrival (TOAs) used in 
Splaver et al. 2005 and those of the later observations, for a 
data span of 21 years, with a noticeable gap between 1994 
and 1998, during the Arecibo upgrade. We used daily av¬ 
eraged TOAs from Splaver et al. 2005, because the original 
190 s integration TOAs are not accessible. Data timestamps 
are derived from observatory masers and retrocorrected to 
Universal Coordinated Time (UTC) via GPS and then fur¬ 
ther corrected to the TT(BIPM) timescale using the 2012 ver¬ 
sion BIPM clock corrections with extrapolations to 2013. The 
TOAs are measured from the observational data through a se¬ 
ries of steps. First the data are folded, as they were being 
taken, into pulse profiles using an ephemeris known to be 
good enough for predicting the pulse period for the duration of 
the observation. The folded profiles from different frequency 
channels and sub-integrations are often summed together to 
improve the signal-to-noise ratio of the profiles. Orthogo¬ 
nal polarizations are summed to produce a total-intensity pro¬ 
file. The summed profiles are then compared with a well- 
measured standard pulse profile from the appropriate fre¬ 
quency band. We employ frequency-domain cross-correlation 
techniques (Taylor 1992) to determine the phase of the pulse 
peak relative to the midpoint of the observation. The final 
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Table 1 

21 year J1713+0747 observations. 


System 

Alias* 

Observatory 

Dates 

Number of 
ToAs 

Epochs 

Bandwidth 

(MHz) 

Typical Integration 
Time (minutes) 

Mark III-At( 1410 MHz) 

M3A-L 

Arecibo 

1992 Jun-1993 Jan 

9 

9 

40 

47 

Mark HUB* (1410 MHz) 

M3B-L 

Arecibo 

1993 Jan-1994 Jan 

46 

46 

40 

47 

Mark IV (1410 MHz) 

M4-L 

Arecibo 

1998 Jul—2004 May 

81 

81 

10 

58 

Mark IV (2380 MHz) 

M4-S 

Arecibo 

1999 Oct-2004 May 

44 

44 

10 

29 

Mark IV-O* (1410 MHz) 

M40-L 

Arecibo 

2004 Jun-2005 Mar 

22 

16 

10 

60 

Mark IV-O* (2380 MHz) 

M40-S 

Arecibo 

2004 Jun—2005 Jan 

8 

7 

10 

30 

ABPP (1410 MHz) 

ABPP-L 

Arecibo 

1998 Feb-2004 May 

98 

89 

56 

60 

ABPP (2380 MHz) 

ABPP-S 

Arecibo 

1999 Dec-2004 May 

46 

46 

112 

30 

ASP (1410 MHz) 

ASP-L 

Arecibo 

2005 Jan—2012 Jan 

990 

48 

64 

20 

ASP (2350 MHz) 

ASP-S 

Arecibo 

2005 Jan-2012 Mar 

668 

41 

64 

20 

GASP (800 MHz) 

GASP-8 

GBT 

2006 Mar-2011 Jan 

997 

41 

64 

20 

GASP (1410 MHz) 

GASP-L 

GBT 

2006 Mar-2010 Jun 

863 

42 

64 

20 

GUPPI (800 MHz) 

GUPPI-8 

GBT 

2010 Mar-2013 Oct 

3533 

49 

800 

20 

GUPPI (1410 MHz) 

GUPPI-L 

GBT 

2010 Mar-2013 Nov 

4381 

64 

800 

20 

PUPPI (1410 MHz) 

PUPPI-L 

Arecibo 

2012 Mar-2013 Nov 

1972 

26 

800 

20 

PUPPI (2300 MHz) 

PUPPI-S 

Arecibo 

2012 Mar-2013 Nov 

992 

24 

800 

20 


These short names are used in Figure 1 and 2 and Table 4. 
f Filter bank used a 78 //.s time constant, 
t Filter bank used a 20 //s time constant. 

* Flere Mark IV-O stands for the recently processed Mark IV data that partially overlap with ASP data. 


TOA of a summed profile is then calculated by adding the 
mid-observation time and the product of pulse period and the 
measured peak phase. The flux density of the pulsar in these 
observations can also be measured by comparing the signal 
strength in the data with that of a calibration observation taken 
right before or after the pulsar observation in which a sig¬ 
nal with known strength was injected. For the post-upgrade 
Arecibo data and all the GBT data, the flux density of the cal¬ 
ibration signal is calibrated every month by comparing it with 
an astronomical object of known and constant flux density, in 
this case, the AGNs J1413+1509 and B1442+09. 

2.3. Instrumental offsets 

The telescopes and data acquisition systems introduce vary¬ 
ing degrees of computational and electronic delays into the 
measured TOAs. Further, in many cases, different standard 
pulse profiles were used to create TOAs from different data 
acquisition systems. As a result, the TOAs from different sys¬ 
tems differ by small time offsets. When a transition is made 
from one system to another at a telescope, typically data are 
collected in parallel during a period of around a year, and 
those data are used to measure the offset between data taken 
with the instruments. Here we describe the measured offsets 
in more detail. 

There was no overlap between data taken with Mark III 
(through 1994) and data taken with Mark IV and ABPP (be¬ 
ginning in 1998). The offset between these two systems was 
treated as a free parameter when fitting timing solutions to the 
full data set. 

Mark IV and ABPP were used in parallel for J1713+0747 
observations between 1998 and 2004. They collected data 
with different bandwidths (Table 1) and were computed us¬ 
ing different profile templates. To align the ABPP ToAs with 
Mark IV, we fitted a phase offset between the 1410-MHz 
TOAs of the two instruments, and found that ABPP ToAs trail 
those of Mark IV by 0.46791 ±0.00009 in pulse phase. In the 
timing modeling, we fix the phase offsets between Mark IV 
and ABPP to this value in both 1410 MHz and 2300 MHz, 
and fit an extra time offset for 2300 MHz ABPP TOAs to ac¬ 
count for any extra template misalignment. 


Mark IV/ABPP and ASP were used in parallel for several 
epochs. We fit across the overlap data to determine the offset 
between the 1410 MHz TOAs of these data sets to measure 
an offset of 2.33 ±0.10 //s (Mark IV trailing ASP). The offset 
between the 2300 MHz TOAs of Mark IV and ASP is treated 
as a free parameter in the full timing solution. 

For the offsets between ASP and PUPPI at Arecibo, and 
GASP and GUPPI at the GBT, analysis of simultaneous ob¬ 
servations of many pulsars, including both pulsar signals and 
radiometer noise, were used to measure very precise offsets 
between the instruments at each observatory; details will be 
given in a forthcoming paper (Arzoumanian et al. 2015). Fur¬ 
ther, the same standard pulse profiles were used in any given 
band. Thus these data sets form a continuous 9-year collec¬ 
tion of TOAs with no arbitrary offsets. 

Finally, the offset between the Arecibo 1410 MHz TOAs 
and GBT Bank 1410 MHz TOAs was treated as a free pa¬ 
rameter when fitting timing solutions to the full data set. The 
offset between the 800 MHz and 1410 MHz GBT TOAs and 
the offset between the 1410 MHz and the 2300 MHz Arecibo 
TOAs are also fitted as free parameters. 

The date span, number of observation epochs, and specifi¬ 
cations of the systems are listed in Table 1. 

3. TIMING MODEL 

We employed the pulsar timing packages TEMPO 17 and 
TEMP02 (Hobbs et al. 2006) to model the TOAs. 

The rotation of the pulsar was modeled with a low-order 
polynomials expansion of spin frequency v and v (Tables 2 
and 3) in order to account for the pulsar’s spin and spin down. 

The pulsar’s position ( a , S ) and proper motion (ji a , fif) on 
the sky and its parallax 03 were also measured through tim¬ 
ing modeling. The distance inferred from our parallax mea¬ 
surement is Dpsr = 1.15 ±0.03 kpc. This distance is con¬ 
sistent with the VLBA parallax distance of 1.05 ±0.06 kpc 
(Chatterjee et al. 2009). 

We employed the Damour & Deruelle (1986) (DD) model 
of binary motion to fit for the binary parameters, including the 

17 http : / /tempo . source forge . net 
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mass of the white dwarf (M c ) measured through Shapiro delay 
(Section 4.1), the orbital period F\,, angle of periastron u>, time 
of periastron passage 7o, projected semi-major axis x and its 
change rate x due to proper motion of the orbit. We observed 
an apparent x as the projection angle of the orbit changed over 
time due to the perpendicular part of the binary’s motion to 
our line of sight. This allowed us to determine the orienta¬ 
tion of the orbit in the sky when combined with the system’s 
proper motion. The orientation of the orbit in the sky is mod¬ 
eled by the parameter f l, the position angle of the ascending 
node. In TEMPO, we grid search for the best f l and then 
hold it fixed when fitting other parameters (Table 2). In the 
T2 model of TEMP02 (Edwards et al. 2006), f l is explicitly 
modeled and fitted, while the changing of x, including the x 
caused by proper motion and the periodic changes due to or¬ 
bital parallaxes of the Earth and the pulsar (Kopeikin 1996), 
are implicitly modeled and not fitted as a parameter (Table 3). 

In our model, we fixed the orbit’s periastron advance rate 
oj to the value inferred from GR and the best-fit binary pa¬ 
rameters (Table 2, 3). This was done by iteratively updating 
the values of oj and refitting for new binary parameters many 
times until the results converged. Fitting for w as a free pa¬ 
rameter resulted in a best-fit value > 3a away from the GR 
prediction. This is likely because, in J1713+0747’s nearly 
circular orbit, <Jj is strongly covariant with the orbital period 
Pd- 

Compared with previous timing efforts, we detected, for the 
first time, an apparent change in the binary period P^, which 
we attribute to the motion of the binary system relative to the 
Sun. This is described in Section 4.2. 

The DMX model was used to fit dispersion measure (DM) 
variations caused by changes in the ISM along the line of sight 
(see Section 3.2 for details). The FD model was used to model 
profile evolution in frequency (see Section 3.4 for details). 

In order to account for unknown systematics in TOAs from 
different instruments, and observation-correlated noise such 
as pulse jitter noise from pulsar emission process, we em¬ 
ployed a general noise model that parameterizes both uncor¬ 
related and correlated noise. The noise modeling is discussed 
in Section 3.1 and the noise model parameters are listed in 
Table 4. 

We used the JPL DE421 solar system ephemeris 
(Folkner et al. 2009) to remove pulse time-of-flight variation 
within the solar system. This ephemeris is oriented to the 
International Celestial Reference Frame (ICRF) and thus our 
astrometric results are also given in the ICRF frame. We 
note parenthetically that a previous generation ephemeris, 
DE405, gave nearly identical but marginally better timing fits 
(Ax 2 - 6 for 14528 dof in both TEMPO and TEMP02). 

The timing parameters and uncertainties are calculated us¬ 
ing a generalized least square (GLS) approach available in 
the TEMPO and TEMP02 software packages. Furthermore 
we also run a Markov Chain Monte Carlo (MCMC) (simi¬ 
lar to Lentati et al. 2014; Vigeland & Vallisneri 2014) anal¬ 
ysis that simultaneously includes the noise parameters (See 
Section 3.1) and the nonlinear timing model. As shown in Ta¬ 
ble 3 and Figure 5, the results from this analysis are very con¬ 
sistent with the GLS approach, indicating that the assumption 
of linearity holds over the full timing parameter range. 

Using a new noise modeling technique, we detected a sig¬ 
nificant red noise signal that could be the same “timing noise” 
described in Splaver et al. (2005). We also detected for the 
first time a change in the observed orbit period. The new tim¬ 
ing model parameters (Table 2 and 3) changed slightly from 


those in Splaver et al. (2005) but are consistent within their 
reported uncertainties. 

3.1. Noise model 

The noise model used in this analysis is a parameter¬ 
ized model that is a function of several unknown quan¬ 
tities describing both correlated and uncorrelated noise 
sources. Uncorrelated noise is independent from one TOA 
to another, while the correlated noise is not. For in¬ 
stance, the template matching error mostly due to radiome¬ 
ter noise are uncorrelated in time, but the pulse jitter 
noise (Cordes & Shannon 2010), which affects the multi¬ 
frequency TOAs measured simultaneously across the band, 
is correlated in time. There is also time-correlated noise, 
such as red timing noise that is correlated from epoch to 
epoch. Among the various types of noise only the tem¬ 
plate matching error a can be estimated when we compute 
the TOAs. Other sources of noise must be modeled sep¬ 
arately. The solution to this problem has been discussed 
extensively in van Haasteren & Levin (2013); Ellis (2013); 
van Haasteren & Vallisneri (2015, 2014); Arzoumanian et al. 
(2014) and Ellis (2014). In this section, we summarize the 
noise model used in this analysis (see e.g. Arzoumanian et al. 
2015, for more details). 

To model uncorrelated noise, we use the standard EFAC and 
EQUAD parameters for each backend/receiver system (e.g. 
PUPPI backend with L-band receiver). These parameters sim¬ 
ply re-scale the original TOA uncertainties 

° 7 ,* —> Ek{a^ k + Qf) ^ , ( 1 ) 

where E\,_ and Q/, denote the EFAC and EQUAD parameters, 
respectively, and the subscript i is the TOA number and the 
subscript k denotes the backend/receiver system. 

To model correlated noise we use the new ECORR param¬ 
eter and a power-law red noise spectrum. The ECORR term 
describes short timescale noise that is completely correlated 
for all TOAs in a given observation but completely uncor¬ 
related between observations. This term could be described 
as pulse phase jitter (Cordes & Shannon 2010) but could also 
have other components. The ECORR term manifests itself as 
a block diagonal term in the noise covariance matrix where 
the size of the blocks is equal to the number of TOAs in a 
given observation. The exact details of the implementation are 
described in Arzoumanian et al. (2015); however, the term es¬ 
sentially acts as a observation-to-observation variance. Lastly, 
we model the red noise as a stationary gaussian process that 
is parameterized by a power spectrum of the form 

( r \ 7red 

UJ ’ (2) 

where A le d is the amplitude of the red noise in //s yr 1 / 2 , 7 rc(l is 
the spectral index. 

These noise parameters are included in a joint likelihood 
that contains all timing model parameters. For the purposes 
of noise modeling, we analytically marginalize over the linear 
timing model parameters and explore the space of noise pa¬ 
rameters via MCMC. We then use the MCMC results to de¬ 
termine the maximum likelihood noise parameters which are 
subsequently used as inputs to TEMPO/TEMP02 GLS fit¬ 
ting routines. In our noise model we include EFAC, EQUAD, 
and ECORR parameters for data collected by different back¬ 
end and receiver systems. However, we do not model EFAC 
and ECORR of the Mark III, Mark IV, and ABPP data because 
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Table 2 

Timing model parameters 3 from TEMPO. 


Parameter 

EFAC and EQUAD 

With Jitter Model 

Jitter and Red Noise Model 

Measured Parameters 

R.A., a (J2000) 

17:13:49.5320251(5) 

17:13:49.5320248(7) 

17:13:49.5320252(8) 

Decl., S (J2000) 

7:47:37.506131(12) 

7:47:37.506155(19) 

7:47:37.50614(2) 

Spin frequecy is (s -1 ) 

218.81184385472585(6) 

218.81184385472594(10) 

218.8118438547251(9) 

Spin down rate is (s -2 ) 

-4.083889(4) x 10“ 16 

-4.083894(7) x 10~ 16 

-4.08382(5) x 10~ 16 

Proper motion in a, = ol cos 5 (mas yr -1 ) 

4.9177(11) 

4.9179(18) 

4.917(2) 

Proper motion in 5, = 5 (mas yr -1 ) 

-3.917(2) 

-3.915(3) 

-3.913(4) 

Parallax, 03 (mas) 

0.858(15) 

0.84(3) 

0.85(3) 

Dispersion measure* 1 (pc cm -3 ) 

15.9700 

15.9700 

15.9700 

Orbital period, P^ (day) 

67.82513682426(16) 

67.82513826935(19) 

67.82513826930(19) 

Change rate of P\,. Pb (10 12 s s~*) 

0.23(12) 

0.41(16) 

0.44(17) 

Eccentricity, e 

0.0000749394(3) 

0.0000749399(6) 

0.0000749402(6) 

Time of periastron passage, Tq (MJD) 

53761.03227(11) 

53761.0328(3) 

53761.0327(3) 

Angle of periastron 0 , lj (deg) 

176.1941(6) 

176.1967(15) 

176.1963(16) 

Projected semi-major axis, x (lt-s) 

32.34242243(5) 

32.34242188(14) 

32.34242188(14) 

sin i, where i is the orbital inclination angle 

0.9672(11) 

0.951(4) 

0.951(4) 

Companion mass, M c ( Mq ) 

0.233(4) 

0.287(13) 

0.286(13) 

Apparent change rate of x, jc (lt-s s -1 ) 

0.00637(7) 

0.00640(10) 

0.00645(11) 

Profile frequency dependency parameter, FD1 

-0.00016317(19) 

-0.0001623(2) 

-0.00016(3) 

Profile frequency dependency parameter, FD2 

0.0001357(3) 

0.0001350(3) 

0.00014(3) 

Profile frequency dependency parameter, FD3 

-0.0000664(6) 

-0.0000668(6) 

-0.000067(17) 

Profile frequency dependency parameter, FD4 

0.0000147(4) 

0.0000153(4) 

0.000015(5) 

Fixed Parameters 

Solar system ephemeris 

DE421 

DE421 

DE421 

Reference epoch for a, 6, and is (MJD) 

53729 

53729 

53729 

Solar wind electron density no (cm _3 ) 

0 

0 

0 

Rate of periastron advance, uj (deg yr -1 ) d 

0.00020 

0.00024 

0.00024 

Position angle of ascending node, Ct (deg) e 

88.43 

88.43 

88.43 

Red noise amplitude (/ is yr 1 / 2 ) 

- 

- 

0.025 f 

Red noise spectral index, 7 re d 

- 

- 

-2.92 

Derived Parameters 

Intrinsic period derivative, Pi nt (s s -1 )* 

8.966(12) x 10“ 21 

8.98(2) x 10“ 21 

8.97(2) x 10- 21 

Dipole magnetic field, B (G) 

2.0485(14) x 10 s 

2.050(3) x 10 s 

2.049(3) x 10 s 

Characteristic age, r c (year) 

8.076(11) x 10 9 

8.07(2) x 10 9 

8.07(2) x 10 9 

Pulsar mass, Mpsr (Mq) 

0.97(3) 

1.32(11) 

1.31(11) 


a We used a modified DD binary model (Damour & Deruelle 1986) that allows us to assume a position angle of ascending node (Q) and fit for the 
apparent change rate of the projected semi-major axis (x) due to proper motion. Numbers in parentheses indicate the 1 cr uncertainties on the last digit(s). 
Uncertainties on parameters are estimated by the TEMPO program using infoimation in the covariance matrix, 
k The averaged DM value; see Section 3.2 and Figure 2 for more discussion. 
c See Figure 2 of Splaver et al. 2005 for definition. 

^ The rate of periastron advance was not fitted but fixed to the GR value because it is highly co-variant with the orbital period. 
e We optimized 12 using a grid search and held it fix to the value that minimized x - 
' The value corresponds to 8.7 X 10~ 15 in the dimensionless strain amplitude unit. 

These parameters are corrected for Shklovskii effect and Galactic differential accelerations. 


there were not enough TOAs in the legacy data set to constrain 
both EFAC and EQUAD. Furthermore, there was only one 
TOA per observation so we cannot constrain the observation- 
correlated noise modeled by ECORR. Instead, we set EFAC 
values to 1 and ECORR to 0 for these data sets, and use only 
EQUAD to model the white noise in these observations. The 
maximum likelihood values of the white noise parameters are 
presented in Table 4. 

Shannon & Cordes (2012) studied the pulse arrival times 
from a single long exposure of PSR J1713+0747, and found 
that this pulsar’s single pulses showed random jitter of ~ 
26 [is. A similar result of ~ 27 /rs was found by Dolch et al. 
(2014) from a more recent study using a 24 hr continuous 
observation of PSR J1713+0747 conducted with major tele¬ 
scopes around the globe. Therefore, by averaging many 
pulses collected in the typical ~ 20 min utes NANOGrav ob¬ 
servation, one expects ~ 27 fis / \/l200v = 51 ns of jitter noise. 
Tables 2 and 3 show the best-fit timing parameters before and 
after we applied our noise model to the data. It is clear that 
the jitter-like observation-correlated noise affected the arrival 


time of the pulses, such that some timing parameters changed 
significantly after including the jitter model. The optimal 
jitter parameters (ECORR, as shown in Table 4) from our 
noise modeling are mostly consistent with the prediction from 
Shannon & Cordes (2012), with some of them being higher. 
This could be due to the covariance between the jitter param¬ 
eters and the EQUAD parameters. 

In Figure 1 we show the red noise realization based on our 
best noise model (Table 4) and compare it to the post-fit resid¬ 
uals of a TEMPO GLS fit. The bottom panels of Figure 1 
show the one- and two-dimensional posterior probability plots 
of the red noise. This noise model describes the data well as 
we can see in Figure 2 in which the maximum likelihood re¬ 
alizations of red and jitter noise are subtracted out. We see 
from the figure that both the high and low frequency residuals 
(with red and jitter noise realizations subtracted) are white 
(described by our EFAC and EQUAD parameters) and the 
weighted residuals follow a zero-mean unit-variance Gaus¬ 
sian distribution. We do note that the normalized residuals do 
not seem to follow exactly the gaussian distribution outside of 
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Table 3 

Timing model parameters 3 from TEMP02. 


Parameter 

EFAC & EQUAD 

With Jitter Model 

Jitter & Red Noise Model 

Measured Parameters 

R. A., a (J2000) 

17:13:49.5320254(5) 

17:13:49.5320247(7) 

17:13:49.5320261(10) 

Decl., 5 (J2000) 

7:47:37.506130(13) 

7:47:37.506130(19) 

7:47:37.50615(3) 

Spin frequecy v (s -1 ) 

218.81184385472573(6) 

218.81184385472589(10) 

218.8118438547255(5) 

Spin down rate v (s -2 ) 

-4.083883(5) x 10~ 16 

-4.083892(7) x 10~ 16 

-4.08386(5) x 10~ 16 

Proper motion in a, = d; cos 5 (mas yr _I ) 

4.9161(12) 

4.9181(18) 

4.915(3) 

Proper motion in <5, fis = 8 (mas yr -1 ) 

-3.915(2) 

-3.910(3) 

-3.914(5) 

Parallax, G5 (mas) 

0.872(16) 

0.83(3) 

0.87(3) 

Dispersion measure* 3 (pc cm -3 ) 

15.9700 

15.9700 

15.9700 

Orbital period, Pj, (day) c 

67.8251365449(12) 

67.8251383194(16) 

67.8251383185(17) 

Change rate of P\,, Pj, (1CT 12 s s -1 ) 

0.19(13) 

0.39(16) 

0.36(17) 

Eccentricity, e 

0.0000749395(3) 

0.0000749400(6) 

0.0000749402(6) 

Time of periastron passage, Tq (MJD) 

53761.03208(11) 

53761.0327(3) 

53761.0328(3) 

Angle of periastron, lj (deg) 

176.1930(6) 

176.1963(15) 

176.1966(14) 

Projected semi-major axis, x (lt-s) 

32.34242258(5) 

32.34242189(13) 

32.34242187(13) 

Orbital inclination, i (deg) 

76.1(3) 

71.9(7) 

71.9(7) 

Companion mass, M c (Mq) 

0.222(4) 

0.286(13) 

0.286(12) 

Position angle of ascending node, O (deg) 

74.3(14) 

89.6(20) 

88(2) 

Profile frequency dependency parameter, FD1 

-0.0001634(2) 

-0.0001628(2) 

-0.0001628(2) 

Profile frequency dependency parameter, FD2 

0.0001358(3) 

0.0001355(3) 

0.0001355(3) 

Profile frequency dependency parameter, FD3 

-0.0000658(6) 

-0.0000671(6) 

-0.0000672(6) 

Profile frequency dependency parameter, FD4 

0.0000141(4) 

0.0000154(4) 

0.0000155(4) 

Fixed Parameters 

Solar system ephemeris 

DE421 

DE421 

DE421 

Reference epoch for a, S, and v (MJD) 

53729 

53729 

53729 

Solar wind electron density uq (cm _3 ) 

0 

0 

0 

Rate of periastron advance, to (deg yr _1 ) d 

0.00019 

0.00024 

0.00024 

Red noise amplitude (/xs yr 1 / 2 ) 

- 

- 

0.025 e 

Red noise spectral index, 7 re d 

- 

- 

-2.92 

Derived Parameters 

Intrinsic period derivative, Pi n t(s s -1 )* 

8.957(13) x 10- 21 

8.98(2) x 10~ 21 

8.96(2) x 10~ 21 

Dipole magnetic field, B (G)* 

2.0473(15) x 10 8 

2.050(3) x 10 s 

2.048(3) x 10 s 

Characteristic age, r c (year)* 

8.085(12) x 10 9 

8.06(2) x 10 9 

8.08(2) x 10 9 

Pulsar mass, Mpsr (Mq) 

0.90(3) 

1.31(11) 

1.31(11) 


a We used TEMP02’s T2 binary model, which implicitly account for the changes of the projected semi-major axis, including x due to proper motion of the 
binary (this allows us to fit for the position angle of ascending node, Q) and the changes due to the orbital parallaxes of the earth and the pulsar (Kopeikin 
1996; Edwards et al. 2006). Numbers in parentheses indicate the 1 a uncertainties on the last digit(s). Uncertainties on parameters are estimated by the 
TEMP02 program using information in the covariance matrix. These uncertainties are consistent with the MCMC results using the full non-linear timing 
model (Section 3.1), see Figure 5 for examples. 

b The averaged DM value; see Section 3.2 and Figure 2 for more discussion. 

c TEMP02’s T2 model reports the orbital period after correcting for the changes in periastron angle to due to proper motion and orbital parallaxes. 
TEMPO DD model does not account for this and therefore reports a slightly different value. 

d The rate of periastron advance was not fitted but fixed to the GR value because it is highly co-variant with the orbital period. 
e This value corresponds to 8.7 X 10 -15 in the dimensionless strain amplitude unit. 

These parameters are corrected for Shklovskii effect and Galactic differential accelerations. 


the 3-sigma range, this affects only < 0.3% of the TOAs and 
will not significantly affect the results presented here. 

Red noise was previously reported for PSR J1713+0747 by 
Splaver et al. (2005), who modeled it using an eighth-order 
polynomial. The largest residual in the present data set ap¬ 
pears in the time period 1999 to 2005 (Figure 1). This is when 
dual-frequency observations begin, and it is likely that the red 
noise model is absorbing unmodeled DM variations in singe- 
frequency data collected before 1999, thus making physical 
interpretations of red noise difficult. 

We reanalyzed the Splaver et al. 2005 data sets (Mark III, 
Mark IV, and ABPP data) with the new red noise model¬ 
ing technique, and show that the timing results from the new 
method are mostly consistent with those from Splaver et al. 
2005, except for a, S, v and i), which probably changed due 
to the new red noise model, and Keplerian parameters, which 
probably changed due to the fact that we used TEMP02’s T2 
model instead of TEMPO’S DD model used by Splaver et al. 
(2005). We fit the Keplerian and post-Keplerian parameters 


simultaneously using the T2 model, whereas, Splaver et al. 
(2005) fit only for the Keplerian parameters while hold¬ 
ing the post-Keplerian parameters to their best-fit values us¬ 
ing the DD model. Therefore, the uncertainties reported in 
Splaver et al. 2005 do not reflect the covariance between the 
two sets of parameters but ours do. 

The red noise signal found by our noise model applied to 
the Splaver et al. (2005) data set is consistent with that found 
in the 21 year data set, both in terms of the noise parame¬ 
ters (Table 2, 3, and 5), and in terms of the shapes of the 
red noise realization (Figure 1), while the red noise modeled 
by high-order frequency polynomials varies significantly de¬ 
pending on the order of the polynomial or and the observation 
time span. 

3.2. DM variation 

The DM of a pulsar reflects the number of free electrons be¬ 
tween the pulsar and the telescopes and it varies because our 
sight-line through the turbulent ISM and solar wind is chang- 
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Figure 1. Top panel: daily-averaged residuals of J1713+0747 produced from a GLS that includes a full noise model. The red dashed line and the gray shaded 
area show the maximum likelihood red noise realization and one-sigma uncertainty. The thick dashed vertical lines separate out various generations of backends 
used at both AO and GBT. Bottom-left panel: two-dimensional marginalized posterior probability plot of red noise spectral index vs. logarithm of the red noise 
amplitude where the solid, dashed and dashed-dotted lines represent the 50, 90, and 95 percent credible regions. The “x” denotes the maximum likelihood value 
of the spectral index and amplitude. Bottom-right panel: one-dimensional marginalized posterior probability for the red noise spectral index where the solid line 
denotes the maximum marginalized a-posteriori value and the dashed lines denote the 68% credible interval. Note that the ID maximum marginalized posterior 
spectral index usually differs slightly from the global maximum likelihood value (—2.92) due to correlations with other parameters. 


Table 4 

Noise parameters 3 and residual rms in //.s. 


Backends 

CT b 

EFAC 

EQUAD 

ECORR 

WRMS C 

AWRMS d 

All 

0.927 




0.246 

0.092 

M3A-L 

0.267 


0.599 


0.589 

0.588 

M3B L 

0.167 


0.412 


0.434 

0.432 

M4-L 

0.172 


0.153 


0.365 

0.146 

M4-S 

0.183 


0.357 


0.668 

0.416 

M40-L 

0.416 


0.315 


0.324 

0.112 

M40-S 

0.355 


0.008 


0.277 

0.141 

ABPP L 

0.106 


0.154 


0.288 

0.067 

ABPP-S 

0.134 


0.260 


0.464 

0.303 

ASP-L 

0.512 

0.979 

0.035 

0.105 

0.222 

0.073 

ASP-S 

0.631 

1.149 

0.004 

0.127 

0.257 

0.115 

GASP-8 

1.391 

1.178 

0.000 

0.023 

0.589 

0.098 

GASP-L 

0.966 

1.128 

0.040 

0.037 

0.288 

0.080 

GUPPI-8 

1.550 

1.052 

0.086 

0.204 

0.657 

0.156 

GUPPI-L 

0.855 

1.204 

0.025 

0.054 

0.266 

0.046 

PUPPI L 

0.303 

1.160 

0.001 

0.094 

0.145 

0.075 

PUPPI S 

0.653 

1.050 

0.058 

0.114 

0.189 

0.080 


11 The unmodeled EFAC and ECORR values default to 1 and 0, respectively. 
b The averaged TOA template matching errors (o). 

c Here WRMS is defined as where is the tim¬ 

ing residual of TOA i, w/ = 1 /of is the weight determined by the TOA errors 
including EFAC and EQUAD, and & — (^2 ^2 w l is the weighted mean 

of the residuals after removing a red-noise model (as in Figure 1 ). 
d AWRMS stands for the weighted RMS of epoch-averaged residuals. 


ing as the pulsar, the Sun, the Earth, and the ISM all move 
with respect to each other. DM variation can affect the timing 
of high-precision pulsars significantly. 

We fit simultaneously with other parameters the time- 
varying DM using the DMX model in TEMPO. This model 
fits independent DM values for TOA groups taken within 14 
day intervals, except for the L-band-only Mark III TOAs. We 
grouped the Mark III TOAs together as a single group, be¬ 
cause their frequency resolution and timing precision are not 
sufficient for measuring epoch-to-epoch DM changes. 

We turned off the solar wind model for TEMPO and 
TEMP02 by setting the solar wind electron density (at 1 AU 
from the Sun) parameter no to 0 cm -3 (the default value is 
10 cm -3 ), and used the DMX model to model all DM varia¬ 
tions including contribution from solar wind. 

Figure 3 shows the measured DM variation of 
PSR J1713+0747. The sudden dip and recovery of DM 
around 2008 (MJD 54800) is due to changes either in the 
ISM or in the solar wind. This DM dip is also observed 
independently by the Parkes observatory (Keith et al. 2013) 
and the European Pulsar Timing Array (G. Desvignes et. al. 
2015, in preparation). 

Spectrum analysis of the time variation of flux, pulse arrival 
phase, and DM have been employed to study the turbulent 
nature of the ISM (e.g. Cordes et al. 1986; Rickett & Lyne 
1990). It has been shown that DM variations of some pulsars 
are consistent with those expected from an ISM characterized 
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Figure 2. Top panel: low frequency daily averaged residuals of J1713+0747 with maximum likelihood jitter and red noise realization subtracted out. Bottom 
panel: high frequency daily averaged residuals of J1713+0747 with maximum likelihood jitter and red noise realization subtracted out. The insets in both top and 
bottom panels show the histogram of the weighted residuals (weighted by EFAC and EQUAD corrected TOA uncertainties) in logarithmic scale in blue curves, 
along with a zero-mean unit-variance Gaussian distribution marked as red dashed curves. 



Year 


Figure 3. The plot shows 16-years DM variation of PSR J1713+0747. The 
dotted line shows the Solar elongation of the pulsar. The subplot shows the 
structure function (error bars) and its a power law fit (solid line). The best-fit 
power law index is 0.49(5), different from the value of 5/3 expected from a 
“pure” Kolmogorov medium. 


by a Kolmogorov turbulence spectrum (Cordes et al. 1990; 
Rickett 1990; Kaspi et al. 1994; You et al. 2007; Keith et al. 
2013; Fonseca et al. 2014). One can calculate the structure 
function of the varying DM: 

0*(r)= (jf) ([ DM(t + r)-DM(t )] 2 ), (3) 

where r is a given time delay, K = 4.148 x 
10 3 MHz 2 pc 1 cm 3 s, and / is the observing frequency in 
MHz. We expect, under the simplest assumptions, this func¬ 


tion to follow a Kolmogorov power law D^(t) = (t/to)' 3-2 , 
where (3 = 11/3 and to is a characteristic time scale related 
to the inner scale of the turbulence. The pulsars with DM 
variations that fit this theory generally have large DM 
variations on timescale of years. However, PSR J1713+0747 
does not show significant long-term DM variation (Figure 3). 
Conversely, it went through a steep drop and recovery around 
2008. If such rapid DM changes are the result of variations 
in the ISM along light of sight, such ISM variations do not fit 
the general characteristics of a Kolmogorov medium. 

3.3. Pulsar spin irregularity 

The term “timing noise” in pulsar timing generally refers 
to the non-white noise left in the timing residuals. An im¬ 
portant contribution to timing noise is expected to come from 
the pulsar’s spin irregularity, i.e., its long-term deviation from 
a simple linear slow down. Spin irregularity is often signifi¬ 
cant in younger pulsars, and may be modeled with high-order 
frequency polynomials (such as v, where v is the pulsar’s 
spin frequency). Potential causes of irregular spin behavior 
include unresolved micro-glitches, internal superfluid turbu¬ 
lence, magnetosphere variations, or external torques caused 
by matter surrounding the pulsar (Hobbs et al. 2010; Yu et al. 
2013; Melatos & Link 2014). These mechanisms could lead 
to accumulative random perturbations in the pulsar’s pulse 
phase, spin rate, or spin-down rate. Shannon & Cordes (2010) 
pointed out that one could model these types of timing noise 
using random walks. Random walks in phase (RWo) would 
grow over time ( T ) proportionally to 7 ■ 2 , random walks in 
v grow proportionally to T 3 / 2 , random walks in i) grow pro¬ 
portionally to T 5 / 2 . Such spin noise would likely have a steep 







































































































































21 years timing of PSR J1713+0747. 


9 


Table 5 

The Timing Results from Splaver et al. 2005 and from a Re-analysis of the Splaver et al. (2005) 
Data set Using New Red Noise Analysis Technique. 


Parameter Splaver et al. 2005 Red Noise Model a 


R. A., ex. (J2000) 

Decl., 8 (J2000) 

Spin frequecy v (s -1 ) 

Spin down rate v (s~ 2 ) 

Proper motion in a, = dcos6 (mas yr -1 ) 
Proper motion in 8, = 8 (mas yr -1 ) 

Parallax, 05 (mas) 

Dispersion measure (pc cm -3 ) 

Orbital period, P\> (day) 

Change rate of P^, /\> (10 -12 s s _1 ) 
Eccentricity, e 

Time of periastron passage, Tq (MJD) 

Angle of periastron, u (deg) 

Projected semi-major axis, x (lt-s) 

Cosine of inclination, cos i 
Companion mass, M c (Mq) 

Position angle of ascending node, Q (deg) 
Solar system ephemeris 
Reference epoch for a, 8, and v (MJD) 
Pulsar mass, Mpsr (Mq) 

Red noise amplitude (/is yr 1 / 2 ) 

Red noise spectral index 


17:13:49.5305335(6) 

17:13:49.5305321(6) 

7:47:37.52636(2) 

7:47:37.52626(2) 

218.8118439157321(3) 

218.811843915731(1) 

-4.0835(2) x 10~ 16 

-4.0836(1) x 10~ 16 

4.917(4) 

4.917(4) 

-3.93(1) 

-3.93(1) 

0.89(8) 

0.84(4) 

15.9960 

15.9940 

67.8251298718(5) b 

67.825129921(4) 

0.0(6) 

-0.2(7) 

0.000074940(1 ) b 

0.000074940(1) 

51997.5784(2) b 

51997.5790(6) 

176.192(l) b 

176.195(3) 

32.34242099(2) b 

32.3424218(3) 

0.31(3) 

0.32(2) 

0.28(3) 

0.30(3) 

87(6) 

89(4) 

DE405 

DE405 

52000 

52000 

1.3(2) 

1.4(2) 

- 

0.004 

- 

5.14 


a We used TEMP02’s T2 binary model, which models the Keplerian (A,, x, e. To, and uj) and post-Keplerian 
orbital elements (cos /', f2, and m 2 ) simultaneously. 

k Splaver et al. 2005 uses TEMPO’S DD model and reports the uncertainties of the Keplerian parameters with 
the post-Keplerian ones fixed to their bestfit values. 


power spectrum with more power in the lower frequencies. 
This is considered as one of the main sources of “red” noise 
in pulsar timing. 

The timing noise of radio pulsars has been studied 
by Cordes & Helfand (1980); Cordes & Downs (1985); 
Arzoumanian et al. (1994); D’Alessandro et al. (1995); 
Matsakis et al. (1997), and later by Hobbs et al. (2010) and 
Shannon & Cordes (2010) with large samples. Matsakis et al. 
(1997) adopted a generalized Allen Variance (traditionally 
used in measuring clock stability) to characterize the timing 
instability of pulsars: 

^(r)=^=(c 2 ) 1/2 , (4) 

where (c 2 ) denotes the sum of squares of the cubic terms fitted 
to segments of length r. Hobbs et al. (2010) found a best- 
fit scaling model of <7,(10 yr) from a large sample of pulsars, 
including canonical pulsars (CPs) and MSPs: 

log 10 M10yr)] = -1.371og 10 [t/ a2 V| a55 ] + 0.52, (5) 

where u, v are the pulsar’s spin and spin-down rate. We find 
that Hobbs et al. (2010)s scaling model (cr™^ 6 * ~ 1 x 10" 12 ) 
over-predicted cr” e ^™ ed = 5 x 10“ 16 for PSR J1713+0747 by 
more than three orders of magnitude. 

Cordes & Helfand (1980) defined a different timing noise 
characteristic <t 2 n 2 based on the root mean square of residuals 
2 f rom a timing fit that does not include any higher order 

spin parameters like v. The timing noise term is related to 

2 

a ^,2 : 

a Tn,2^T) + <J wi ( 6 ) 

where rr^, is a time-independent term caused by white noise 
in the data. In this definition, timing noise cr 2 N2 (r) grows 


bigger over time while white noise stays constant. 

Shannon & Cordes (2010) studied the (Jj N2 from a large 
sample of CPs and MSPs. They found a scaling model: 

ln(a XN , 2 )= 1.6—1.41n(i/)+l.lln|t>_ 15 |+21n(r yr ), (7) 

where zx _ 15 is i> in units of 10 -15 s -2 , and T yl is the observation 
time span in years. This scaling model predicts that, for 21- 
year timing of PSR J1713+0747, the residual RMS without 
removing timing noise cr£ N2 would be ~ 400 ns. The mea¬ 
sured RMS of the red noise residual a 1 ^ RN = 364 ns, is con¬ 
sistent with the extrapolation from Shannon & Cordes (2010). 
The best-fit scaling law also indicates that the residuals of the 
sampled pulsars <jtn ,2 seem to grow linearly with T yr . If the 
timing noise of the sampled pulsars is due to the accumulation 
of spin noise, and the spin noise is caused by the same phys¬ 
ical processes, then this RMS growth rate would imply that 
the spin noise of pulsars has a frequency power spectrum of 
power-law index 7 rc . (i ~ -5. This spectral index is consistent 
with the 7 re d from our noise model. This can be seen in the 
bottom right plot of Figure 1 by noting that a spectral index 
of -5 is consistent with the posterior at the one-sigma level. 

It is inconclusive whether or not the observed red noise can 
be interpreted as pulsar spin irregularity. Other sources of 
noise also could have contributed significantly. If we do as¬ 
sume that they are from spin irregularity, the estimated max¬ 
imum likelihood red noise spectral index of ~ -3 favors that 
the pulsar spin irregularities come from random walks in ei¬ 
ther spin phase or spin rate, although other explanations can¬ 
not be ruled out due to the substantial uncertainty on the red 
noise spectral index (Figure 1). 

Finally, Shannon & Cordes (2010) showed that the signifi¬ 
cance of timing noise coming from gravitational wave (GW) 
background could be estimated as crew ,2 ~ 1.34o(T yr ) 5/3 ns. 
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Figure 4. GUPPI and PUPPI post-fit residuals versus frequency when fitted 
without the FD model, showing the frequency dependence of the TOAs that 
is not accounted for by the DM. For clarity, this plot is made using only TOAs 
with error smaller than 3 fis. 

where Aq is the characteristic strain at / = 1 year -1 and T yr is 
the observational time span in years. The current best upper 
limit on GW characteristic strain is 2.4 x 10~ 15 (Shannon et al. 
2013), which predicts an upper limit on timing noise of ~ 
500 ns from GW background. Therefore, we cannot rule out 
the contribution of gravitational waves in the timing noise. 

3.4. Pulse profile evolution with frequency 

After removing the dispersion that causes TOA delays 
proportional to f~ 2 , where / is the observing frequency, 
we still see small remaining frequency-dependent residuals 
from wide-band observations using different instruments and 
telescopes (Figure 4). It appears that the low-frequency 
(~800 MHz) signals lead the high-frequency (1400 and 
2300 MHz band) signals by microseconds. The cause of such 
TOA evolution is not clear. It could be a change in the pul¬ 
sar’s radiation pattern with frequency, or it could be the use 
of different standard profiles in different frequencies. Pulsar 
radiation of different frequencies may originate from different 
parts of the star’s magnetosphere, and the radiation region of 
the pulsars’ magnetosphere may be slightly distorted, leading 
to a frequency-dependent radiation pattern. Pennucci et al. 
(2014) and Liu et al. (2014) extensively discussed this phe¬ 
nomenon and developed TOA extraction techniques based on 
phase-frequency 2-D pulse profile matching. This technique 
is not yet applied to our data set. 

Demorest et al. (2013) allowed an arbitrary offset between 
TOAs taken with different observing systems and at differ¬ 
ent frequencies in order to model profile evolution with fre¬ 
quency. However, the number of frequency channels has in¬ 
creased by a factor of ten with the modern wide-band instru¬ 
ments, making it much harder to mitigate profile-frequency 
evolution using frequency channel offsets. Instead, we used 
the FD model, a polynomial of the logarithm of frequency: 
AfpD = XX i c,(log/)' (Arzoumanian et al. 2015; solid line in 
Figure 4) to fit for and remove the profile-frequency evolu¬ 
tion, where AfpD is the profile evolution term in unit of sec¬ 
ond, / is the observing frequency in unit of GHz and c, are 
the FD model parameters. We employed an F-test with sig¬ 
nificance value of 0.0027 to determine how many FD param¬ 
eters are needed to model profile frequency evolution. PSR 
J1713+0747 only requires n = 4 FD parameters. 

4. RESULTS 


4.1. Mass measurements 

The timing model of PSR J1713+0747 has been signifi¬ 
cantly improved by the 21-year timing effort. Most notably, 
the pulsar and the companion masses have been more pre¬ 
cisely constrained (Table 2 and 3, Figure 5) through Shapiro 
delay measurements. The companion’s mass M c = 0.286 ± 
0.012 Mq and the pulsar Mpsr = 1.31 ±0.11 M© are in good 
agreement with the previously measured values (Splaver et al. 
2005). Furthermore, we have carried out an MCMC run that 
uses the nonlinear timing model in order to map out any 
non-Gaussian correlations in parameter space. We find that 
the nonlinear model gives nearly identical results to the GLS 
method of TEMPO/TEMP02. The covariance matrix used 
in our GLS fitting contains terms come from both the corre¬ 
lated and the uncorrelated noise; therefore, the timing param¬ 
eter uncertainties we get have taken into account the contri¬ 
bution from the noise model. We note that, without the noise 
model, the derived pulsar masses would be substantially, and 
perhaps unrealistically, lower (with AMpsr ranges from 0.3 
to 0.4 Mq) than the values with the noise model (Table 2 and 
3), suggesting that correlated noise would significantly impact 
the accuracy of high precision timing analysis. 

The pulsar’s mass is compatible with the distribution of pul¬ 
sar masses in other neutron star-white dwarf systems, and 
in good agreement with the distribution of pulsar masses 
found in recycled binaries (Ozel et al. 2012; Kiziltan et al. 
2013). The precise measurement of neutron star masses 
may eventually help us understand the properties of mat¬ 
ter of extreme density (Demorest et al. 2010; Lattimer 2012; 
Antoniadis et al. 2013). 

In the standard picture of binary evolution, an MSP with 
a low-mass white dwarf companion must have been spun up 
through accretion when the white dwarf was a giant star fill¬ 
ing its Roche lobe. This should lead to a strong correla¬ 
tion between the binary period and the mass of the white 
dwarf companion (Rappaport et al. 1995; Tauris & Savonije 
1999; Podsiadlowski et al. 2002). Indeed, this picture has 
been supported the measurements of several pulsar bi¬ 
nary systems (e.g. van Straten et al. 2001; Kaspi et al. 1994; 
Tauris & van den Heuvel 2014; Ransom et al. 2014). The or¬ 
bital period and companion mass of PSR J1713+0747 fit this 
correlation very well, thus supporting the standard MSP evo¬ 
lution theory. 

4.2. Intrinsic orbital decay 

We have observed an apparent change in orbital period from 
PSR J1713+0747, P h = (0.36 ±0.17) x 10“ 12 s s -1 (Tables 
2 and 3). This orbital period change is not intrinsic to the 
pulsar binary, but rather the result of the relative accelera¬ 
tion between the binary and the observer, i.e. the combination 
of differential acceleration in the Galactic gravitational po¬ 
tential (Damour & Taylor 1991) and the “Shklovskii” effect 
(Shklovskii 1970) which is caused by the transverse motion 
of the pulsar binary relative to Earth. We have good measure¬ 
ments of the distance and proper motion of the binary system, 
which allow us to remove these effects and study the system’s 
intrinsic orbital decay. 

The apparent change in orbital period due to differential ac¬ 
celeration in Galactic gravitational potential can be derived 
from 

pGat _ = (-0.10 ±0.02) x 10 -12 s s -1 , (8) 

c 

where Aq is the line of sight acceleration of the pulsar binary. 
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Figure 5. One and two-dimensional posterior probability distributions of the cosine of the inclination angle, pulsar mass, and companion mass from the noise 
model MCMC including the full nonlinear timing model. The maximum marginalized posterior value and \a credible interval is very consistent with the GLS 
solution from TEMPO/TEMP02. The solid, dashed and dashed-dotted lines represent the one, two and three-sigma credible regions, respectively. 


Ag is obtained using Equation 5 in Nice & Taylor (1995), 
Equation 17 in Lazaridis et al. (2009), the local matter den¬ 
sity of Galactic disk around solar system (Holmberg & Flynn 
2004) and the Galactic potential model by Reid et al. (2014). 

The Shklovskii effect causes to change by 

pf k = (^ 2 +/ 4)-P b = (0.65 ±0.02) x 10“ 12 ss _1 . (9) 

c 

Therefore, the pulsar’s intrinsic orbital decay is = P° bs - 
pShk_pGai _ (-0.20±0.17) x 10 -12 s s -1 , and is consistent 
with zero. 

Due to the very long ~68 day orbit, the binary’s decay due 
to the emission of gravitational radiation is expected to be un¬ 
detectable: P// R = -6 x 10 -18 s s -1 (Lorimer & Kramer 2005). 
Therefore, the insignificant intrinsic orbital decay rate is en¬ 
tirely consistent with the description of quadrupolar gravita¬ 
tional radiation within General Relativity (GR). 

Other than the gravitational radiation, two classical effects 
could have played a role in P b nt . One, /© is caused by mass 
loss in the binary system, and the other, p£, is the contribution 
from tidal effects. The pulsar and the white dwarf both could 
lose mass due to their magnetic dipole radiation; the maxi¬ 
mum mass loss rate due to this effect can be estimated from 
the star’s rotational energy loss rate. In the case of the pulsar, 
Mpsr = E /c 2 , measurable through the spin down rate of the 
pulsar. The white dwarf generally loses mass at a much lower 
rate than the pulsar. Therefore, orbital change due to mass loss 
can be estimated as P/f ~ 1 x 1 (T 14 s s -1 (Damour & Taylor 
1991; Equation (9) and (10) of Freire et al. 2012b). This is an 
order of magnitude smaller than the measured uncertainties 
on The tidal effect in this binary system is expected to 
be P b -C 3 x 10‘ 14 s s -1 based on the most extreme scenarios 
(the white dwarf spins at its break-up velocity and the tidal 
synchronizing time scale equals the characteristic age of the 
pulsar; see Equation 11 in Freire et al. 2012b and references 


therein). Both of these extra terms are much smaller than the 
observed uncertainties on . 

4.3. Time Variation of G 

Based on the measurement of “excess” orbital pe¬ 
riod change Pf c = Pjf - Pf - P/ - P® R , Damour et al. 
(1988) derived a phenomenological limit for G without 
considering the binding energy of the stars: G/G ~ 
_pe xc /( 2 p b ) _ (1.0 ±2.3) x 10 -11 yr -1 using the timing of 
binary PSR B1913+16. Since then P“ c of pulsar binaries, 
including PSR J1713+0747, have been used to constrain 
G/G (Kaspi et al. 1994; Nice et al. 2005; Deller et al. 2008; 
Lazaridis et al. 2009; Freire et al. 2012b). So far all pulsar ob¬ 
servations show G/G consistent with being zero, with upper 
limits largely determined by the uncertainties in orbital period 
change rate, distance, and proper motions. PSR J1713+0747 
has the smallest known P“ c /(2Pb) — (-0.5 ±0.9) x 10 -12 yr -1 
(Section 4.2) and is particularly useful for constraining the 
time variability of the gravitational constant. 

Nordtvedt (1990), Lazaridis et al. (2009), and Freire et al. 
(2012b) showed that a generic test of G/G can be achieved 
using pulsar binaries in a more rigorous fashion by incorpo¬ 
rating the binding energy of the neutron stars. The binding 
energy of a compact star changes with the G, resulting in a 
changing mass, this will also affect the binary orbit. In a 
generic form, we could characterize this effect using a self¬ 
gravity "sensitivity" parameter s p (Will 1993). The changing 
G will now change the orbital period of a pulsar binary system 
(Nordtvedt 1990; Lazaridis et al. 2009): 



This formalism is more generic in the sense that it incorpo¬ 
rates the compactness of the neutron star, but it also assumes 
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Figure 6. Confidence contour of G/G and kb calculated from 
PSRs J1012+5307, J1738+0333, and J1713+0747 using an MCMC sim¬ 
ulation. The shaded area marks the 95% confidence G limit from LLR 
(Hofmann et al. 2010). The gray area marks the 95% confidence G limit 
from planetary ephemerides (Fienga et al. 2014). 


that non-perturbative effects are absent and higher order con¬ 
tributions in the self-gravity sensitivity can be neglected. 

Meanwhile, in the framework of an alternative gravitation 
theory that violates SEP, a binary system may emit dipole 
gravitational radiation (Will 1993, 2001; Lazaridis et al. 2009; 
Freire et al. 2012b and references therein). Such effects arise 
when the two bodies are very different in terms of their self¬ 
gravity, i.e. their compactness. Under the aforementioned as¬ 
sumptions of neglecting non-perturbative effects and higher 
order contributions of self-gravity sensitivity, this extra dipole 
radiation could lead to an extra orbital change term: 

P? - -471^ k d S 2 , (11) 

(Will 1993; Lazaridis et al. 2009), where Tq = GM 0 /c 3 = 
4.925490947 /is, fj, is the reduced mass ( m p m c /M ) of the sys¬ 
tem, «'/) is a dipole gravitational radiation “coupling constant,” 
and S is the difference between the self-gravity “sensitivities” 
of the two bodies (S = s p -s c \ s p ~ 0.1 m p /Mo according to 
Damour & Esposito-Farese 1992 ; and s c <C s p ). In Einstein’s 
GR kd = 0 — there is no self-gravity induced dipole gravita¬ 
tional radiation, but it is generally not the case in alternative 
theories that violate the SEP. 

PSR J1713+0747 has a wider binary orbit than most other 
high-timing-precision pulsar binaries, making its P® very 
small. Conversely, is larger when is large. This makes 
PSR J1713+0747 the best pulsar binary system for constrain¬ 
ing the effect of the changing gravitational constant G. Lim¬ 
its on both G and Kd can be estimated in the same fashion 
as in Lazaridis et al. (2009): by solving G and kd simulta¬ 
neously from the equation P“ c = P// + P// (Equation 29 of 
Lazaridis et al. 2009) of different pulsars. We applied this 
method to four pulsars: PSR J0437-4715, PSR J1012+5307, 
PSR J1738+0333, and PSR J1713+0747 using timing param¬ 
eters reported in Lazaridis et al. (2009), Freire et al. (2012b), 
and this work. The resulting confidence region of G and kd 


is shown in Figure 6. We found, at 95% confidence limit, 
G/G = (Obi 1.1) x 1(T 12 yr 1 ; k d = (-0.9±3.3) x KT 4 
This constraint on G is more stringent than previous pulsar- 
based constraints (Freire et al. 2012b), and close to one of 
the best constraints of this type ( G/G = (-0.07 ± 0.76) x 
10 -12 yr -1 ) from the Lunar Laser Ranging (LLR) experi¬ 
ment (Hofmann et al. 2010), which measured Earth-Moon 
distance to ~ 10“’ 1 precision through 39 years of laser rang¬ 
ing. Fienga et al. (2014) showed that G/G can be constrained 
to (0.01±0.18)x 10 _12 yr _1 through the analysis of solar sys¬ 
tem planetary ephemerides. The Kd limit, which is not con¬ 
strained by solar-system tests, is also slightly improved by 
using PSR J1713+0747. The pulsar-timing G and Kd limits 
are particularly interesting in the testing of SEP-violating al¬ 
ternative theories, because they arise from a test using objects 
of strong self-gravitation. For example, in some classes of the 
scalar-tensor theories, the effect of G/G could be significantly 
enhanced in pulsar-white dwarf binaries (Wex 2014). 


4.4. SEP and PFE 


The SEP states that the gravitational effect on a small test 
body is independent of its constitution, and in particular, that 
bodies of different self-gravitation should behave the same 
in the same gravitational experiments. This principle is vi¬ 
olated in alternative theories of gravitation like the aforemen¬ 
tioned Jordan-Fierz-Brans-Dicke scalar-tensor theory. The 
PSR J1713+0747 binary is an excellent laboratory for test¬ 
ing effects of SEP violation. If the SEP is violated, the neu¬ 
tron star and the white dwarf will be accelerated differently 
by the Galactic gravitational field, causing the binary orbit to 
be polarized toward the center of the Galaxy. The excess ec¬ 
centricity is expected to be (Damour & Schafer 1991): 


k+l = 7T- 


A gic 2 


- (T21 

2-F£(MpsR + M c )(27r/ J P b ) 2 ’ 
where J 7 is a factor accounts for potential changes in the peri- 
astron advance rate due to deviations from GR, Q is the effec¬ 
tive gravitational constant in the interaction between the pul¬ 
sar and the white dwarf, and g± is the projection of Galactic 
acceleration on the orbital plane and A is the dimensionless 
factor that characterizes the significance of SEP violation. We 
assume TQ ~ G here and after. The Galactic acceleration of 
the pulsar system is derived from Holmberg & Flynn (2004); 
Reid et al. (2014). 

GR predicts that there is no preferred reference frame in the 
universe but this may be different in alternative theories. The 
Parameterized Post-Newtonian (PPN) formalism parameter¬ 
ized possible deviations from GR into a set of parameters (see 
Will 2014 for the list of them), some of which are associated 
with the PFE. In this work, we test 0:3, one of the PFE-related 
parameters. If a 3 B 0, this would lead to both the presence of 
a PFE and the breaking of momentum conservation. A rotat¬ 
ing body would be accelerated perpendicular to its spin axis 
and its absolute velocity in the preferred reference frame. In a 
pulsar binary, this effect would cause an excess in eccentricity, 
which can be estimated by (Damour & Esposito-Farese 1992; 
Bell & Damour 1996): 


= «3 


c p\ w \Pb 


-sin/3, 


(13) 


24nP TQ(M PS r+M c ) 
where w is the absolute velocity of the binary system rela¬ 
tive to the preferred frame of reference, typically taken as 
that of the cosmic microwave background (CMB), P is the 
pulsar’s spin period, B is the angle between w and the spin 
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axis of the pulsar, and A3 is the strong-field version of the 
PPN parameters <23. Here w = w Q + v P sr, where w Q = 369 ± 
0.9 km s _1 is the velocity of the solar system relative to the 
CMB (Hinshaw et al. 2009), and the term vpsr is the relative 
speed of the pulsar to our solar system. Vpsr is only partially 
known because we can measure the pulsar system’s proper 
motion on the sky but we cannot measure its line of sight ve¬ 
locity (l'r). 

Fortunately, many variables in these equations are measur¬ 
able in the case of the PSR J1713+0747 binary. It is possible 
to constrain A and A3 using Bayesian techniques by assuming 
certain fiducial priors for v r (Splaver et al. 2005; Stairs et al. 
2005; Gonzalez et al. 201 1). In our case, we assumed a Gaus¬ 
sian prior for v r centred at zero with a width equal to the sys¬ 
tem’s proper motion speed. Based on our 21-year timing of 
J1713+0747 alone, we find 95% confidence limits on the vi¬ 
olations of SEP and Lorentz invariance A < 0.01 and A 3 < 
2 x 10 -20 , slightly improving the single pulsar limits from ear¬ 
lier data on this pulsar (Splaver et al. 2005; Stairs et al. 2005; 
Gonzalez et al. 201 1). Stronger limits can be found by com¬ 
bining the results from multiple similar pulsar systems (Wex 
2000; Stairs et al. 2005; Gonzalez et al. 201 1). 

5. SUMMARY 

In this paper, we present a comprehensive model of high 
precision timing observations of PSR J1713+0747 that spans 
21 years. We improved measurements of the pulsar and its 
companion’s masses and the shape and orientation of the bi¬ 
nary orbit. We also detect, for the first time, an apparent 
change in orbital period due to Galactic differential acceler¬ 
ations and the Shklovskii effect. These measurements, when 
combined with those of other pulsars, significantly improve 
the pulsar timing limit on the rate of change of the gravita¬ 
tional constant, G. Although the pulsar constraint is not bet¬ 
ter than the best solar system ones, it is nevertheless an in¬ 
dependent test using extra-solar binary systems thousands of 
light-years away. The pulsar tests also could be more con¬ 
straining for some classes of alternative gravitational theories 
that predict stronger non-GR effects in strong-field regime 
(Wex 2014). The new best pulsar timing limit on G/G is 
(0.6 ± 1.1) x 10 -12 yr -1 (< 0.033 H 0 based on the 3 -a limit), 
where Ho is the Hubble constant. In other words, the change 
rate of gravitational constant has to be a factor of at least 31 
(3-cr limit) slower than the average expansion rate of the Uni¬ 
verse. 

Meanwhile, the precise measurements of PSR 
J1713+0747’s orbital eccentricity and 3D orientation al¬ 
low us to test the violation of SEP and Lorentz invariance 
with it. We found a single-pulsar 95% upper limit on 
A < 0.01, the SEP violation factor, and A3 < 2 x 10 -20 , 
the PPN parameter that characterizes violation of Lorentz 
invariance. Because of the different statistical analysis 
methods used, our A and A3 limits are slightly different but 
still consistent with the results of the same tests in previous 
publications (Wex 2000; Splaver et al. 2005; Stairs et al. 
2005; Gonzalez et al. 2011; Freire et al. 2012a). Ultimately, 
the newly discovered pulsar triple system PSR J0337+1715 
(Ransom et al. 2014) could yield the best test on SEP viola¬ 
tion (Freire et al. 2012a; Shao et al. 2015; Berti et al. 2015). 
In this case the inner pulsar-white dwarf binary is orbited by 
another white dwarf in an outer orbit, making this system an 
excellent laboratory for testing the free fall of a neutron star 
and white dwarf in an external gravity field. 

We studied the time variation of PSR J1713+0747’s DM 


from 1998 to 2013, and fitted the structure function of the 
DM variation with a power law. The best-fit power law in¬ 
dex is 0.49(5), significantly smaller than the 5/3 index ex¬ 
pected from a “pure” Kolmogorov medium. This relatively 
flat structure function could be the result of either the lack 
of long-term DM variations or an excess of short-term vari¬ 
ations. The sudden DM dip around 2008 (Figure 3) is a 
good example of such short-term DM variations. Similar non- 
Kolmogorov DM variations have been observed in some of 
the other NANOGrav pulsars (L. Levin et al. 2015, in prepa¬ 
ration.). Evidence of non-Kolmogorov behavior in the ISM 
was also found in the analysis of multi-frequency pulsar scat¬ 
ter times (Lewandowski et al. 2015). 

As part of our timing modeling, we also included noise con¬ 
tribution such as jitter and red noise using the GLS fit and a 
covariance matrix that included the correlated and uncorre¬ 
lated noise terms. We found that our timing result is signifi¬ 
cantly affected by the noise model, especially the jitter noise, 
suggesting that the adoption of jitter modeling may be neces¬ 
sary in other cases of high precision pulsar timing. We found 
that our noise parameters and timing residuals are consistent 
with the jitter noise estimates from Shannon & Cordes (2012) 
and the timing noise estimate from Shannon & Cordes (2010). 
However, the scaling law extrapolated from large sample stud¬ 
ies of timing noise in Hobbs et al. (2010) overestimated the 
timing noise level a : ( 1 Oyr) in this pulsar. 

Our noise model parameters and timing residual RMS (Ta¬ 
ble 4) provide a crude estimation of the amount of noise in 
our data. The weighted root mean square (WRMS; see Table 
4 for definition) of the 21-year daily-averaged timing residu¬ 
als is ~ 92 ns. Table 4 shows a systematic improvement in 
the timing accuracy of this pulsar in the last two decades, due 
to advances in instrumentation. But the improvements are not 
as large as expected from the radiometer equation, perhaps 
because of pulse jitter. 

Assuming that the red noise is caused by spin irregularity, 
the best-fit spectral index is consistent with spin irregularity 
caused by random walks in either the spin phase or the spin 
rate of the pulsar, but it does not exclude other explanations 
due to its large uncertainty (see bottom panel of Figure 1). 
The observed red noise level is also consistent with the pre¬ 
diction from the current best upper limit of GW background 
(Shannon et al. 2013), therefore, we cannot rule out signifi¬ 
cant timing noise contribution from the GW background. 
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